The Dynamics of Supercooled Silica: Acoustic modes and Boson 

peak 

t^. | Jiirgen Horbach, Walter Kob and Kurt Binder 

Institut fur Physik, Johannes Gutenberg- Universitdt, Staudinger Weg 7, D-55099 Mainz, 

Germany 
(February 1, 2008) 

■ Using molecular dynamics computer simulations we investigate the dy- 

namics of supercooled silica in the frequency range 0.5-20 THz and the wave- 
vector range 0.13-l.lA -1 . We find that for small wave-vectors the dispersion 
r~| ! relations are in very good agreement with the ones found in experiments and 

that the frequency at which the boson-peak is observed shows a maximum at 
g ■ around 0.39A" 1 . 
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PACS numbers: 61.43.Fs, 61.20.Lc, 02.70.Ns, 64.70.Pf 



I. INTRODUCTION 

In the last few years a significant effort was undertaken to understand the nature of the 
so-called boson-peak, a prominent dynamical feature at around 1 THz which is observed 
in strong glassformers JT] ^] . Various theoretical approaches have been proposed to explain 
this peak, such as localized vibrational modes or scattering of acoustic modes, but so far no 
clear picture has emerged yet. Recently also computer simulations have been used in order 
■ to gain insight into the mechanism that gives rise to this peak, but due to the high cooling 
rates with which the samples were prepared (on the order of 10 12 K/s) and small system sizes 
(20-40A) the results of these investigations were not able to give a final answer either @,f|. 
In the present work we present the results of a large scale computer simulation of supercooled 
silica. At the temperature investigated we are still able to fully equilibrate the sample, thus 
avoiding the problem of the high cooling rates, and by using a large system we can minimize 
the possibility of finite size effects ||. Thus, by making a large computational effort, we are 
able to study the dynamics of this strong glassformer in a frequency and wave-vector range 
which is not accessible to real experiments and can therefore investigate the properties of 
q ■ the boson-peak in greater detail than was possible so far. 
O 

• Eh . II. DETAILS OF THE SIMULATIONS 

X 

The silica model we use for our simulation is the one proposed by van Beest et al. || 
and has been shown to give a good description of the static properties of silica glass as 
well as the dynamical properties of the supercooled melt, such as the activation energy of 
the diffusion constant |§ . In this model the potential between ions i and j is given by 

^;) = ^ + A;e-^-^ . (1) 
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The values of the parameters Aij, and Cy can be found in the original publication ||. 
The non-Coulombic part of the potential was truncated and shifted at 5.5A. The simulations 
were done at constant volume and the density of the system was fixed to 2.3 g/cm 3 . The 
system size was 8016 ions, giving a size of the box of (48.37A) 3 , and the equations of motion 
were integrated over 4 • 10 6 times steps of 1.6fs, thus over a time span of 6.4ns. This time 
is sufficiently long to fully equilibrate the system at 2900K, the temperature considered in 
this study ||. More details on the simulations can be found in Ref. ||. 



III. RESULTS 

In the present work we study the dynamics of the system by means of J^(g, z/) and 
Jr{q, v), the longitudinal and transverse current-current correlation functions for wave- 



vector q at frequency v |T0|. These are defined as the longitudinal and transverse part 



of the current-current correlation function, i.e. 

/oo 
exp(zq ■ [r k (t) - r ; (0)])) , (2) 
-°° ki 

where Ufc(t) is equal to q • Ykif)/q for a = L and equal to q x r k (t)/q for a = T. Note 
that J^(g, z/) is also equal to b ,2 S(q,b , )/q 2 , where S(q, v) is the dynamical structure factor 
as measured in scattering experiments. In the following we will focus on the silicon-silicon 
correlation only, but we have found that the oxygen-oxygen correlation function behave very 
similarly. 

In Fig. 1 we show the frequency dependence of Jl((?, v) for wave- vectors between 0.13A -1 , 
the smallest wave-vector compatible with our box, and 0.8A -1 , a wave-vector which is still 
significantly smaller than the location of the first sharp diffraction peak in S(q), which is 
around 1.6A -1 . From the figure we recognize that this correlation function has a peak at a 
frequency Vh{q) which increases with increasing q and which correspond to the longitudinal 
acoustic modes. A similar picture is obtained for Jr(q, v) @- The fact that also this corre- 
lation function shows an acoustic mode shows that even at this relatively high temperatures 
the system is able to sustain transverse acoustic modes and thus is visco-elastic. 

The boson-peak is seen best in the dynamic structure factor S(q, u), which is shown in 
Fig. 2 for small values of q. From this figure we see that VBp{q), the location of the boson- 
peak, is g-dependent in that it moves from small frequencies for small values of q (dashed 
lines) to a maximum frequency for q « 0.39A -1 (bold dotted line) and then back to small 
frequencies for large values of q (solid lines) (see also inset of Fig. 3). Also included is the 
location of the boson-peak as determined from neutron scattering experiments at 1673K 
which is around 1.5 THz ||. We will comment more on this work below. 

From the wave-vector dependence of Ul and ut we obtain the dispersion relation which 
is presented in Fig. 3. We see that, as expected, for small wave-vectors v L depends linearly 
on q. Also included in the figure is a line with slope cz,=6370m/s, the experimental value of 
the longitudinal sound velocity of silica at around 1600K M. We see that the data points 
for ul for small q are very close to this line and thus we conclude that the sound velocity of 
this system is very close to the one of real silica, thus giving further support for the validity 
of the model potential. 
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For the transverse acoustic modes the agreement between the experiment and the simu- 
lation data is a bit inferior, but still good. We also note that for wave-vectors larger than 
1.4A -1 , a bit less than the location of the first sharp diffraction peak, vi and v? do not 
increase anymore. The reason for this is likely the fact that at this q value the system has 
a quasi-Brillouin zone [[|. 

Also included in the figure is vbp, the location of the boson-peak. We find that for large 
wave-vectors vbp is around 1.8 THz, a value that is a bit larger than the experimental value 
of 1.5 THz reported by Wischnewski et al. at 1673K ||. However, these authors also found 
that vbp increases with increasing temperature and a rough extrapolation of their data for 
vbp t° T=2900K shows that a value of 1.8 THz is quite reasonable, thus giving further 
support for the validity of our model. 

In the inset of Fig. 3 we show the dispersion curves at small values of q. Interestingly we 
observe that VBp(q) shows a maximum at around q ~ 0.39A -1 . This maximum might be 
due to the fact that the mechanism leading to the boson-peak is particularly effective at this 
wave- vector or that there are two mechanisms giving rise to the boson-peak, one dominating 
at small wave- vectors and the second one dominating at larger wave- vectors and that in the 
vicinity of q ~ 0.39A -1 the sum of the contribution of the two mechanisms is largest. 

Finally we mention that at small wave- vectors the curve for VBp(q) seems to join smoothly 
the one for v^. Within the accuracy of our data it is not clear, whether the boson-peak and 
the longitudinal acoustic mode become identical or whether the boson-peak ceases to exist 
for wave-vectors smaller than approximately 0.2A -1 and thus we cannot use this feature to 
exclude one of the possible mechanisms proposed for the boson-peak. 

To summarize we can say that our simulation allows to investigate the dynamics of 
supercooled silica in a wave-vector and frequency range which is not accessible to real ex- 
periments. We find that the dispersion relations for the longitudinal and transverse acoustic 
modes agree very well with the experimental values and that vbp shows a maximum at 
around 0.39A -1 . This q corresponds to a length scale on the order of ISA. Thus we have 
evidence that the mechanism leading to the boson-peak is very effective on the length scale 
of several tetrahedra. 
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FIG. 1. Frequency dependence of the longitudinal current-current correlation for different 
wave-vectors q. 
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FIG. 2. Frequency dependence of the dynamic structure factor for different wave-vectors q. The 
vertical line is the location of the boson-peak at 1673K as determined from the neutron scattering 
experiments [3]. 



5 




FIG. 3. Wave-vector dependence of vl (open circles), ut (filled circles) and ubp (open triangles). 
The bold solid lines are the dispersion relations for the longitudinal and transverse acoustic modes 
(Ref. [3]). Inset: enlargement of the curves at small q. 
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